1. Setup & Recap

This tutorial extends the 1-pool Bayesian tutorial to a 3-pool series model. If you haven’t gone through the 1-pool tutorial yet, start there – it introduces Bayesian inference, likelihood functions, priors, and MCMC from scratch.

In the 1-pool tutorial, we learned that:

  • Two parameters (decomposition rate k and carbon input) control the model
  • Two data points (bulk C stock and bulk \(\Delta^{14}\)C at 2022) constrain those parameters
  • The posterior was fairly tight because 2 parameters vs. 2 data points is a well-determined problem

Now we move to 6 parameters and 6 data points. The problem is formally just as determined, but the parameter space is much harder to explore and the parameters can trade off against each other in complex ways.

library(SoilR)
library(BayesianTools)
library(ggplot2)
library(dplyr)
library(tidyr)
library(patchwork)

set.seed(42)

load("soilRmcmc.RData")

We reuse the same run_and_summarize helper from the 1-pool tutorial, extended to handle any number of parameters.

# 'prior' can be any BayesianTools prior object (uniform, truncated normal, custom).
# If NULL, a uniform prior is created from prior_lower/prior_upper.
run_and_summarize <- function(likelihood_fn, prior = NULL,
                              prior_lower = NULL, prior_upper = NULL,
                              par_names,
                              iterations = 20000,
                              burn_in = 5000, thin = 10) {

  if (is.null(prior)) {
    prior <- createUniformPrior(lower = prior_lower, upper = prior_upper)
  }

  bt_setup <- createBayesianSetup(
    likelihood = likelihood_fn,
    prior      = prior,
    names      = par_names
  )

  mcmc_out <- runMCMC(
    bayesianSetup = bt_setup,
    sampler       = "DREAMzs",
    settings      = list(iterations = iterations, message = FALSE)
  )

  raw_chains <- getSample(mcmc_out, start = 1, coda = TRUE, parametersOnly = FALSE)
  chain_df <- do.call(rbind, lapply(seq_along(raw_chains), function(i) {
    ch <- as.data.frame(as.matrix(raw_chains[[i]]))
    n_par <- length(par_names)
    colnames(ch)[1:n_par] <- par_names
    ch$log_posterior <- ch[, n_par + 1]
    ch <- ch[, c(par_names, "log_posterior")]
    ch$iteration <- seq_len(nrow(ch))
    ch$chain <- factor(i)
    ch
  }))

  samples <- getSample(mcmc_out, start = burn_in, thin = thin)
  colnames(samples) <- par_names
  posterior_df <- as.data.frame(samples)

  map_est <- MAP(mcmc_out, start = burn_in)$parametersMAP
  names(map_est) <- par_names

  summary_tbl <- posterior_df %>%
    pivot_longer(everything(), names_to = "parameter", values_to = "value") %>%
    group_by(parameter) %>%
    summarise(
      median   = median(value),
      lower_90 = quantile(value, 0.05),
      upper_90 = quantile(value, 0.95),
      .groups  = "drop"
    )

  list(
    mcmc_out     = mcmc_out,
    chain_df     = chain_df,
    posterior_df = posterior_df,
    map          = map_est,
    summary      = summary_tbl,
    burn_in      = burn_in,
    thin         = thin
  )
}

2. Why a 3-Pool Model?

The CiPEHR soil profile has three distinct layers: surface organic, deep organic, and mineral. These layers have very different properties – the surface is fresh, carbon-rich litter that decomposes quickly, while the mineral layer is old, stabilized carbon that turns over on centennial timescales.

When we aggregated everything into a single pool in the 1-pool tutorial, we lost this structure. The bulk \(\Delta^{14}\)C is a weighted average that hides how the bomb spike propagates differently through each layer.

Let’s look at the per-layer data.

dat_layers <- datComb %>%
  filter(year == 2009 | treatment == "Control") %>%
  mutate(layer_label = factor(layer,
    levels = c("so", "do", "min"),
    labels = c("Surface organic", "Deep organic", "Mineral")))

knitr::kable(
  dat_layers %>%
    select(year, layer_label, C_stock, C_stock_sd, C14, C_14_sd) %>%
    arrange(year, layer_label),
  digits = 1,
  col.names = c("Year", "Layer", "C stock (kg/m2)", "C stock SD",
                "Delta14C", "Delta14C SD"),
  caption = "Per-layer soil observations from CiPEHR"
)
Per-layer soil observations from CiPEHR
Year Layer C stock (kg/m2) C stock SD Delta14C Delta14C SD
2009 Surface organic 1.9 0.5 151.0 56.5
2009 Deep organic 10.8 2.6 -59.8 129.0
2009 Mineral 14.3 7.8 -240.0 63.9
2022 Surface organic 1.9 0.5 136.2 62.7
2022 Deep organic 10.8 2.4 -20.1 22.5
2022 Mineral 11.8 7.9 -255.8 89.2 
p_layer_c14 <- ggplot(dat_layers, aes(x = year, y = C14,
                                       color = layer_label, shape = layer_label)) +
  geom_point(size = 4) +
  geom_errorbar(aes(ymin = C14 - C_14_sd, ymax = C14 + C_14_sd), width = 0.3) +
  geom_line(linewidth = 0.8) +
  labs(title = expression(paste(Delta^14, "C by Soil Layer")),
       x = "Year", y = expression(paste(Delta^14, "C (per mil)")),
       color = NULL, shape = NULL) +
  scale_color_manual(values = c("coral", "steelblue", "forestgreen")) +
  theme_minimal(base_size = 14) +
  theme(legend.position = "top")

p_layer_cstock <- ggplot(dat_layers, aes(x = year, y = C_stock,
                                          color = layer_label, shape = layer_label)) +
  geom_point(size = 4) +
  geom_errorbar(aes(ymin = C_stock - C_stock_sd, ymax = C_stock + C_stock_sd),
                width = 0.3) +
  geom_line(linewidth = 0.8) +
  labs(title = "Carbon Stock by Soil Layer",
       x = "Year", y = expression(paste("C stock (kg m"^-2, ")")),
       color = NULL, shape = NULL) +
  scale_color_manual(values = c("coral", "steelblue", "forestgreen")) +
  theme_minimal(base_size = 14) +
  theme(legend.position = "top")

p_layer_c14 + p_layer_cstock

Notice how the layers differ:

  • Surface organic has the highest \(\Delta^{14}\)C – it has incorporated the most bomb \(^{14}\)C because it cycles fastest.
  • Mineral has the lowest \(\Delta^{14}\)C – it cycles so slowly that the bomb spike has barely arrived.
  • Deep organic falls in between.

A 1-pool model can’t capture these differences. The 3-pool model treats each layer as a separate compartment, connected in series:

\[\text{Atmosphere} \xrightarrow{\text{input}} \text{Surface organic} \xrightarrow{a_{21}} \text{Deep organic} \xrightarrow{a_{32}} \text{Mineral}\]

3. The 3-Pool Series Model

Parameters

The 3-pool series model has 6 parameters:

Parameter Description Units
k1 Decomposition rate, surface organic yr\(^{-1}\)
k2 Decomposition rate, deep organic yr\(^{-1}\)
k3 Decomposition rate, mineral yr\(^{-1}\)
input Carbon input rate to surface pool gC/m\(^2\)/yr
a21 Transfer rate, surface \(\to\) deep yr\(^{-1}\)
a32 Transfer rate, deep \(\to\) mineral yr\(^{-1}\)

The compartmental matrix

The model follows the same structure as the 1-pool model – \(dC/dt = I + A \cdot C\) – but now \(C\) is a 3-element vector and \(A\) is a \(3 \times 3\) compartmental matrix:

\[A = \begin{pmatrix} -k_1 & 0 & 0 \\ a_{21} & -k_2 & 0 \\ 0 & a_{32} & -k_3 \end{pmatrix}\]

Each diagonal entry \(-k_i\) is the total loss rate from pool \(i\) (decomposition + transfer out). The off-diagonal entries \(a_{21}\) and \(a_{32}\) route carbon downward. This is a “series” topology – no carbon moves back up.

For this to be physically valid, the transfer rates must be smaller than the decomposition rates: \(a_{21} < k_1\) and \(a_{32} < k_2\). Otherwise, more carbon would leave a pool through transfer than through decomposition, which doesn’t make physical sense.

Forward model demo

Let’s run the model at a few hand-picked parameter sets to build intuition.

dat_3pool <- datComb %>%
  filter(year == 2009 | treatment == "Control") %>%
  mutate(C_stock    = C_stock * 1000,
         C_stock_sd = C_stock_sd * 1000)

init_3pool  <- dat_3pool %>% filter(year == 2009)
final_3pool <- dat_3pool %>% filter(year == 2022)

C0_3pool <- init_3pool$C_stock
F0_3pool <- init_3pool$C14
years_3pool <- seq(2009, 2022)

pool_labels <- c(so = "Surface organic", do = "Deep organic", min = "Mineral")
years_plot <- seq(2009, 2022, by = 0.5)

demo_params <- list(
  list(ks = c(0.1, 0.01, 0.005), input = 100, a21 = 0.05, a32 = 0.005,
       label = "Fast surface"),
  list(ks = c(0.02, 0.01, 0.005), input = 50, a21 = 0.01, a32 = 0.002,
       label = "Slow surface"),
  list(ks = c(0.05, 0.02, 0.01), input = 80, a21 = 0.02, a32 = 0.01,
       label = "High transfer")
)

demo_runs <- bind_rows(lapply(demo_params, function(dp) {
  m <- tryCatch(
    ThreepSeriesModel14(
      t = years_plot,
      ks = setNames(dp$ks, c("k1", "k2", "k3")),
      C0 = C0_3pool, F0_Delta14C = F0_3pool,
      In = dp$input, a21 = dp$a21, a32 = dp$a32,
      inputFc = atmIn, pass = TRUE
    ),
    error = function(e) NULL
  )
  if (is.null(m)) return(NULL)
  data.frame(
    year = rep(years_plot, 3),
    pool = rep(pool_labels, each = length(years_plot)),
    C14     = as.vector(getF14(m)),
    C_stock = as.vector(getC(m)),
    scenario = dp$label
  )
}))

demo_runs$pool <- factor(demo_runs$pool, levels = pool_labels)

obs_plot <- dat_3pool %>%
  mutate(pool = factor(pool_labels[layer], levels = pool_labels))

p_demo_c14 <- ggplot(demo_runs, aes(x = year, y = C14, color = scenario)) +
  geom_line(linewidth = 1) +
  geom_point(data = obs_plot, aes(x = year, y = C14),
             inherit.aes = FALSE, size = 3) +
  geom_errorbar(data = obs_plot,
                aes(x = year, ymin = C14 - C_14_sd, ymax = C14 + C_14_sd),
                inherit.aes = FALSE, width = 0.3) +
  facet_wrap(~ pool, ncol = 1, scales = "free_y") +
  labs(title = expression(paste("Forward model demo: ", Delta^14, "C")),
       x = "Year", y = expression(paste(Delta^14, "C (per mil)")),
       color = NULL) +
  theme_minimal(base_size = 12) +
  theme(legend.position = "top")

p_demo_cstock <- ggplot(demo_runs, aes(x = year, y = C_stock / 1000, color = scenario)) +
  geom_line(linewidth = 1) +
  geom_point(data = obs_plot, aes(x = year, y = C_stock / 1000),
             inherit.aes = FALSE, size = 3) +
  geom_errorbar(data = obs_plot,
                aes(x = year, ymin = (C_stock - C_stock_sd) / 1000,
                    ymax = (C_stock + C_stock_sd) / 1000),
                inherit.aes = FALSE, width = 0.3) +
  facet_wrap(~ pool, ncol = 1, scales = "free_y") +
  labs(title = "Forward model demo: Carbon stock",
       x = "Year", y = expression(paste("C stock (kg m"^-2, ")")),
       color = NULL) +
  theme_minimal(base_size = 12) +
  theme(legend.position = "top")

p_demo_c14 + p_demo_cstock

None of these hand-picked parameter sets perfectly match the observations in all three layers simultaneously. That’s the challenge – and exactly why we need MCMC to explore the 6-dimensional parameter space.

4. The 3-Pool Likelihood

The likelihood has the same Gaussian structure as the 1-pool version, but now we compare predictions to observations in each of the 3 pools separately. That gives us 6 data points: 3 carbon stocks and 3 \(\Delta^{14}\)C values at 2022.

One important difference: we multiply the reported measurement uncertainty by a factor of 0.4. This is an empirical scaling factor from the original analysis that accounts for model structural error – the fact that the model is a simplified representation of reality. Without this scaling, the likelihood would be too forgiving and the posteriors too wide.

likelihood_3pool <- function(pars) {
  k1 <- pars[1]; k2 <- pars[2]; k3 <- pars[3]
  input <- pars[4]; a21 <- pars[5]; a32 <- pars[6]

  model <- tryCatch(
    ThreepSeriesModel14(
      t = years_3pool,
      ks = c(k1 = k1, k2 = k2, k3 = k3),
      C0 = C0_3pool,
      F0_Delta14C = F0_3pool,
      In = input,
      a21 = a21,
      a32 = a32,
      inputFc = atmIn,
      pass = TRUE
    ),
    error = function(e) return(NULL)
  )

  if (is.null(model)) return(-1e10)

  idx <- length(years_3pool)
  pred_C14    <- getF14(model)[idx, ]
  pred_Cstock <- getC(model)[idx, ]

  ll_C14    <- sum(dnorm(pred_C14 - final_3pool$C14, 0,
                         final_3pool$C_14_sd * 0.4, log = TRUE))
  ll_Cstock <- sum(dnorm(pred_Cstock - final_3pool$C_stock, 0,
                         final_3pool$C_stock_sd * 0.4, log = TRUE))

  return(sum(ll_C14, ll_Cstock, na.rm = TRUE))
}

cat("Test likelihood at example parameters:",
    likelihood_3pool(c(0.1, 0.01, 0.01, 200, 0.1, 0.01)))
## Test likelihood at example parameters: -42.98645

Unlike the 1-pool model, we can’t easily visualize a 6-dimensional likelihood surface. We have to trust the MCMC sampler to explore it for us. This is where good priors and sufficient iterations become critical.

5. Priors

Choosing prior distributions

As we discussed in the 1-pool tutorial, the prior encodes what we believe before seeing the data. With 6 parameters, prior choice matters even more – too wide and the MCMC wastes time in impossible regions; too narrow and you might miss the answer.

We have two broad options:

  • Uniform (flat): Equal probability across a range. Simple, but treats all values as equally likely.
  • Truncated normal (Gaussian): Centered on a best guess with a bell-curve shape, clipped to physical bounds. Encodes the idea that some values are more plausible than others.

We’ll start with uniform priors and then compare to truncated normal priors in Section 9.

The key physical constraints for both prior types:

  • Surface organic decomposes fastest (\(k_1\) is largest)
  • Mineral soil decomposes slowest (\(k_3\) is smallest)
  • Transfer rates must be smaller than decomposition rates (\(a_{21} < k_1\), \(a_{32} < k_2\))

Uniform priors

par_names <- c("k1", "k2", "k3", "input", "a21", "a32")

prior_lower <- c(k1 = 0.002,  k2 = 0.0006, k3 = 0.0002,
                 input = 10,   a21 = 0.0025, a32 = 0.0001)
prior_upper <- c(k1 = 0.2,    k2 = 0.06,   k3 = 0.06,
                 input = 205,  a21 = 0.5,    a32 = 0.1)

prior_range_df <- data.frame(
  Parameter = par_names,
  Lower = prior_lower,
  Upper = prior_upper,
  Interpretation = c(
    "Turnover 5-500 yr",
    "Turnover 17-1667 yr",
    "Turnover 17-5000 yr",
    "Litter input to surface pool",
    "Surface to deep transfer",
    "Deep to mineral transfer"
  )
)
knitr::kable(prior_range_df, digits = 4,
             caption = "Uniform prior ranges for the 3-pool model")
Uniform prior ranges for the 3-pool model
Parameter Lower Upper Interpretation
k1 k1 0.0020 0.20 Turnover 5-500 yr
k2 k2 0.0006 0.06 Turnover 17-1667 yr
k3 k3 0.0002 0.06 Turnover 17-5000 yr
input input 10.0000 205.00 Litter input to surface pool
a21 a21 0.0025 0.50 Surface to deep transfer
a32 a32 0.0001 0.10 Deep to mineral transfer

Gaussian (truncated normal) priors

Suppose a literature review suggests approximate values for each parameter. We can encode this as truncated normal priors – Gaussian bell curves centered on the literature values, clipped to the same physical bounds as the uniform priors.

# Literature-informed centers and uncertainties
prior_means <- c(k1 = 0.05,  k2 = 0.015,  k3 = 0.005,
                 input = 80,  a21 = 0.03,   a32 = 0.01)
prior_sds   <- c(k1 = 0.04,  k2 = 0.015,  k3 = 0.01,
                 input = 60,  a21 = 0.05,   a32 = 0.02)

prior_gaussian <- createTruncatedNormalPrior(
  mean  = prior_means,
  sd    = prior_sds,
  lower = prior_lower,
  upper = prior_upper
)

gaussian_range_df <- data.frame(
  Parameter = par_names,
  Mean = prior_means,
  SD = prior_sds,
  Lower = prior_lower,
  Upper = prior_upper
)
knitr::kable(gaussian_range_df, digits = 4,
             caption = "Truncated normal prior parameters (same bounds as uniform)")
Truncated normal prior parameters (same bounds as uniform)
Parameter Mean SD Lower Upper
k1 k1 0.050 0.040 0.0020 0.20
k2 k2 0.015 0.015 0.0006 0.06
k3 k3 0.005 0.010 0.0002 0.06
input input 80.000 60.000 10.0000 205.00
a21 a21 0.030 0.050 0.0025 0.50
a32 a32 0.010 0.020 0.0001 0.10

Comparing prior shapes

# Uniform prior samples
prior_samples_unif <- as.data.frame(mapply(
  runif, n = 10000, min = prior_lower, max = prior_upper
))

# Truncated normal prior samples
prior_samples_gauss <- as.data.frame(mapply(function(mu, s, lo, hi) {
  # Rejection sampling from truncated normal
  out <- numeric(10000)
  filled <- 0
  while (filled < 10000) {
    candidates <- rnorm(20000, mu, s)
    valid <- candidates[candidates >= lo & candidates <= hi]
    n_take <- min(length(valid), 10000 - filled)
    out[(filled + 1):(filled + n_take)] <- valid[1:n_take]
    filled <- filled + n_take
  }
  out
}, mu = prior_means, s = prior_sds, lo = prior_lower, hi = prior_upper))
colnames(prior_samples_gauss) <- par_names

prior_long_unif <- prior_samples_unif %>%
  pivot_longer(everything(), names_to = "parameter", values_to = "value") %>%
  mutate(parameter = factor(parameter, levels = par_names),
         prior_type = "Uniform")

prior_long_gauss <- prior_samples_gauss %>%
  pivot_longer(everything(), names_to = "parameter", values_to = "value") %>%
  mutate(parameter = factor(parameter, levels = par_names),
         prior_type = "Gaussian")

prior_long <- bind_rows(prior_long_unif, prior_long_gauss)

ggplot(prior_long, aes(x = value, fill = prior_type)) +
  geom_density(alpha = 0.4, linewidth = 0.6) +
  facet_wrap(~ parameter, scales = "free", ncol = 3) +
  scale_fill_manual(values = c("Gaussian" = "coral", "Uniform" = "skyblue"),
                    name = "Prior type") +
  labs(title = "Uniform vs. Gaussian (Truncated Normal) Priors",
       subtitle = "Same physical bounds, but Gaussian concentrates probability near the center",
       x = "Parameter value", y = "Density") +
  theme_minimal(base_size = 12) +
  theme(legend.position = "top")

The Gaussian priors have the same physical bounds as the uniform priors, but they concentrate probability near the literature-informed centers. For well-constrained parameters (like \(k_1\)), this may not matter much – the data will drive the posterior regardless. For weakly constrained parameters (like \(k_3\) or \(a_{32}\)), the prior shape can make a real difference.

6. Running MCMC

With 6 parameters, we need more iterations than the 1-pool model. We use 20,000 iterations with DREAMzs. For a teaching example this is sufficient; a publication-quality analysis might use 100,000+.

result <- run_and_summarize(
  likelihood_fn = likelihood_3pool,
  prior_lower   = prior_lower,
  prior_upper   = prior_upper,
  par_names     = par_names,
  iterations    = 20000,
  burn_in       = 5000,
  thin          = 10
)

cat("Retained", nrow(result$posterior_df), "posterior samples")
## Retained 501 posterior samples

Trace plots

With 6 parameters, we need a bigger panel. Look for the same “hairy caterpillar” pattern in each – well-mixed chains bouncing around a stable region.

trace_vars <- c(par_names, "log_posterior")
trace_labels <- c(par_names, "Log-posterior")

chain_long <- result$chain_df %>%
  pivot_longer(cols = all_of(c(par_names, "log_posterior")),
               names_to = "parameter", values_to = "value") %>%
  mutate(parameter = factor(parameter, levels = trace_vars, labels = trace_labels))

ggplot(chain_long, aes(x = iteration, y = value, color = chain)) +
  geom_line(alpha = 0.4, linewidth = 0.3) +
  facet_wrap(~ parameter, scales = "free_y", ncol = 3) +
  labs(title = "MCMC Trace Plots (3-Pool Model)",
       subtitle = "Log-posterior (bottom right) should plateau when chains converge",
       x = "Iteration", y = "Value") +
  theme_minimal(base_size = 11) +
  theme(legend.position = "top")

The log-posterior panel is the most important diagnostic: it shows the objective value being maximized. Look for a clear plateau – the point where the log-posterior stops climbing and starts bouncing around a stable level. This is where the chains have found the high-probability region, and everything before it is burn-in.

Some parameters (like k3 or a32 for the slow mineral pool) may show worse mixing than others. This is typical – slow pools have weak 14C signals, so the data provide less constraint and the chains wander more.

Burn-in and thinning

We use a burn-in of 5,000 (vs. 2,000 for the 1-pool model) because the chains need more time to settle in a 6D space. Thinning of 10 reduces autocorrelation.

ggplot(chain_long, aes(x = iteration, y = value, color = chain)) +
  geom_line(alpha = 0.4, linewidth = 0.3) +
  geom_vline(xintercept = result$burn_in / 3, color = "black",
             linewidth = 1, linetype = "dotted") +
  annotate("text", x = result$burn_in / 3, y = Inf,
           label = " Burn-in cutoff", hjust = 0, vjust = 1.5, size = 3.5) +
  facet_wrap(~ parameter, scales = "free_y", ncol = 3) +
  labs(title = "Burn-in Removal",
       subtitle = "Everything left of the dotted line is discarded. Check that log-posterior has plateaued by the cutoff.",
       x = "Iteration", y = "Value") +
  theme_minimal(base_size = 11) +
  theme(legend.position = "top")

7. Results

Prior vs. posterior

The key question: how much did the data narrow down our estimates for each parameter?

posterior_long <- result$posterior_df %>%
  pivot_longer(everything(), names_to = "parameter", values_to = "value") %>%
  mutate(parameter = factor(parameter, levels = par_names))

ggplot() +
  geom_histogram(data = prior_long_unif,
                 aes(x = value, y = after_stat(density)),
                 bins = 40, fill = "skyblue", alpha = 0.4, color = "white") +
  geom_density(data = posterior_long, aes(x = value),
               fill = "coral", alpha = 0.5, linewidth = 0.8) +
  geom_vline(data = data.frame(
    parameter = factor(par_names, levels = par_names),
    xval = result$map),
    aes(xintercept = xval), linetype = "dashed", color = "darkblue") +
  facet_wrap(~ parameter, scales = "free", ncol = 3) +
  labs(title = "Prior (blue) vs. Posterior (coral)",
       subtitle = "Dashed blue = MAP estimate. Narrow posteriors = well-constrained parameters.",
       x = "Parameter value", y = "Density") +
  theme_minimal(base_size = 12)

MAP estimates and credible intervals

# Compute uncertainty reduction
summary_with_reduction <- result$summary %>%
  mutate(
    prior_range = prior_upper[parameter] - prior_lower[parameter],
    post_range  = upper_90 - lower_90,
    uncert_reduction = round(100 * (1 - post_range / prior_range), 1)
  )

cat("=== 3-Pool Posterior Summary ===\n")
## === 3-Pool Posterior Summary ===
cat(sprintf("%-8s  %8s  %10s  %10s  %10s  %s\n",
            "Param", "MAP", "Median", "90% lo", "90% hi", "Reduction"))
## Param          MAP      Median      90% lo      90% hi  Reduction
for (i in seq_len(nrow(summary_with_reduction))) {
  r <- summary_with_reduction[i, ]
  cat(sprintf("%-8s  %8.4f  %10.4f  %10.4f  %10.4f  %5.1f%%\n",
              r$parameter, result$map[r$parameter], r$median,
              r$lower_90, r$upper_90, r$uncert_reduction))
}
## a21         0.0902      0.1057      0.0695      0.1573   82.4%
## a32         0.0022      0.0119      0.0015      0.0317   69.8%
## input      30.4176     46.9212     16.7631    106.4977   54.0%
## k1          0.0221      0.0323      0.0081      0.0635   72.0%
## k2          0.0165      0.0188      0.0060      0.0348   51.4%
## k3          0.0123      0.0256      0.0040      0.0518   20.1%

Derived quantities: turnover times

The turnover time \(\tau = 1/k\) tells us how long carbon resides in each pool on average. This is often more intuitive than the decomposition rate.

k_params <- c("k1", "k2", "k3")
turnover_df <- result$posterior_df[, k_params] %>%
  mutate(across(everything(), ~ 1 / .x)) %>%
  rename(tau1 = k1, tau2 = k2, tau3 = k3) %>%
  pivot_longer(everything(), names_to = "parameter", values_to = "value") %>%
  mutate(pool = factor(parameter,
    levels = c("tau1", "tau2", "tau3"),
    labels = c("Surface organic", "Deep organic", "Mineral")))

turnover_summary <- turnover_df %>%
  group_by(pool) %>%
  summarise(
    median = median(value),
    lower_90 = quantile(value, 0.05),
    upper_90 = quantile(value, 0.95),
    .groups = "drop"
  )

knitr::kable(turnover_summary, digits = 1,
             col.names = c("Pool", "Median (yr)", "90% lo", "90% hi"),
             caption = "Posterior turnover times (1/k)")
Posterior turnover times (1/k)
Pool Median (yr) 90% lo 90% hi
Surface organic 30.9 15.7 123.2
Deep organic 53.2 28.7 167.4
Mineral 39.1 19.3 248.5

Posterior correlations

In multi-parameter models, parameters often trade off against each other. For example, a higher decomposition rate might be compensated by higher carbon input. These correlations are a key feature of the posterior that single-parameter summaries miss.

cor_mat <- cor(result$posterior_df)

cat("--- Posterior Correlation Matrix ---\n")
## --- Posterior Correlation Matrix ---
print(round(cor_mat, 2))
##          k1   k2   k3 input   a21   a32
## k1     1.00 0.14 0.09  0.84  0.33 -0.29
## k2     0.14 1.00 0.01  0.22  0.42  0.03
## k3     0.09 0.01 1.00  0.09  0.05  0.24
## input  0.84 0.22 0.09  1.00  0.27 -0.23
## a21    0.33 0.42 0.05  0.27  1.00 -0.09
## a32   -0.29 0.03 0.24 -0.23 -0.09  1.00
n_pairs <- min(500, nrow(result$posterior_df))
pairs(result$posterior_df[sample(nrow(result$posterior_df), n_pairs), ],
      pch = 16, cex = 0.4, col = adjustcolor("steelblue", 0.3),
      main = "Posterior Pairwise Correlations (3-Pool Model)")

Look for elliptical clouds tilted at an angle – these indicate correlated parameters. Strong correlations mean the data can’t independently constrain those parameters; only certain combinations are well determined.

8. Model Predictions

We run the forward model at the MAP estimate (best fit) and at 200 random posterior draws (uncertainty envelope), showing predictions for each soil layer separately.

years_plot <- seq(2009, 2022, by = 0.5)

# MAP prediction
map_model <- ThreepSeriesModel14(
  t = years_plot,
  ks = c(k1 = result$map["k1"], k2 = result$map["k2"], k3 = result$map["k3"]),
  C0 = C0_3pool, F0_Delta14C = F0_3pool,
  In = result$map["input"], a21 = result$map["a21"], a32 = result$map["a32"],
  inputFc = atmIn, pass = TRUE
)

map_pred <- data.frame(
  year = rep(years_plot, 3),
  pool = rep(pool_labels, each = length(years_plot)),
  C14     = as.vector(getF14(map_model)),
  C_stock = as.vector(getC(map_model))
)

# Posterior ensemble
n_draws <- min(200, nrow(result$posterior_df))
draw_idx <- sample(nrow(result$posterior_df), n_draws)

ensemble <- bind_rows(lapply(draw_idx, function(i) {
  p <- as.numeric(result$posterior_df[i, ])
  m <- tryCatch(
    ThreepSeriesModel14(
      t = years_plot,
      ks = c(k1 = p[1], k2 = p[2], k3 = p[3]),
      C0 = C0_3pool, F0_Delta14C = F0_3pool,
      In = p[4], a21 = p[5], a32 = p[6],
      inputFc = atmIn, pass = TRUE
    ),
    error = function(e) NULL
  )
  if (is.null(m)) return(NULL)
  data.frame(
    year = rep(years_plot, 3),
    pool = rep(pool_labels, each = length(years_plot)),
    C14     = as.vector(getF14(m)),
    C_stock = as.vector(getC(m))
  )
}))

envelope <- ensemble %>%
  group_by(year, pool) %>%
  summarise(
    C14_lo = quantile(C14, 0.05), C14_hi = quantile(C14, 0.95),
    Cstock_lo = quantile(C_stock, 0.05), Cstock_hi = quantile(C_stock, 0.95),
    .groups = "drop"
  )

obs_plot <- dat_3pool %>%
  mutate(pool = factor(pool_labels[layer], levels = pool_labels))

map_pred$pool  <- factor(map_pred$pool, levels = pool_labels)
envelope$pool  <- factor(envelope$pool, levels = pool_labels)

# Delta14C by layer
p_c14 <- ggplot() +
  geom_ribbon(data = envelope,
              aes(x = year, ymin = C14_lo, ymax = C14_hi),
              fill = "steelblue", alpha = 0.3) +
  geom_line(data = map_pred, aes(x = year, y = C14),
            color = "steelblue", linewidth = 1) +
  geom_point(data = obs_plot, aes(x = year, y = C14), size = 3) +
  geom_errorbar(data = obs_plot,
                aes(x = year, ymin = C14 - C_14_sd, ymax = C14 + C_14_sd),
                width = 0.3) +
  facet_wrap(~ pool, ncol = 1, scales = "free_y") +
  labs(title = expression(paste(Delta^14, "C by Soil Layer")),
       subtitle = "MAP prediction + 90% posterior envelope vs. observations",
       x = "Year", y = expression(paste(Delta^14, "C (per mil)"))) +
  theme_minimal(base_size = 12)

# C stock by layer
p_cstock <- ggplot() +
  geom_ribbon(data = envelope,
              aes(x = year, ymin = Cstock_lo / 1000, ymax = Cstock_hi / 1000),
              fill = "forestgreen", alpha = 0.3) +
  geom_line(data = map_pred, aes(x = year, y = C_stock / 1000),
            color = "forestgreen", linewidth = 1) +
  geom_point(data = obs_plot, aes(x = year, y = C_stock / 1000), size = 3) +
  geom_errorbar(data = obs_plot,
                aes(x = year, ymin = (C_stock - C_stock_sd) / 1000,
                    ymax = (C_stock + C_stock_sd) / 1000),
                width = 0.3) +
  facet_wrap(~ pool, ncol = 1, scales = "free_y") +
  labs(title = "Carbon Stock by Soil Layer",
       subtitle = "MAP prediction + 90% posterior envelope vs. observations",
       x = "Year", y = expression(paste("C stock (kg m"^-2, ")"))) +
  theme_minimal(base_size = 12)

p_c14 + p_cstock

9. Experiment: Uniform vs. Gaussian Priors

Does prior shape matter for the 3-pool model? Let’s re-run the MCMC with the truncated normal priors from Section 5 and compare the posteriors.

result_gauss <- run_and_summarize(
  likelihood_fn = likelihood_3pool,
  prior         = prior_gaussian,
  par_names     = par_names,
  iterations    = 20000,
  burn_in       = 5000,
  thin          = 10
)

cat("Retained", nrow(result_gauss$posterior_df), "posterior samples (Gaussian prior)")
## Retained 501 posterior samples (Gaussian prior)

Comparing posteriors

posterior_long_unif <- result$posterior_df %>%
  pivot_longer(everything(), names_to = "parameter", values_to = "value") %>%
  mutate(parameter = factor(parameter, levels = par_names),
         prior_type = "Uniform")

posterior_long_gauss <- result_gauss$posterior_df %>%
  pivot_longer(everything(), names_to = "parameter", values_to = "value") %>%
  mutate(parameter = factor(parameter, levels = par_names),
         prior_type = "Gaussian")

posteriors_compare <- bind_rows(posterior_long_unif, posterior_long_gauss)

ggplot(posteriors_compare, aes(x = value, fill = prior_type)) +
  geom_density(alpha = 0.4, linewidth = 0.6) +
  facet_wrap(~ parameter, scales = "free", ncol = 3) +
  scale_fill_manual(values = c("Gaussian" = "coral", "Uniform" = "steelblue"),
                    name = "Prior type") +
  labs(title = "Posterior Comparison: Uniform vs. Gaussian Priors",
       subtitle = "Same likelihood and data, different prior shapes",
       x = "Parameter value", y = "Density") +
  theme_minimal(base_size = 12) +
  theme(legend.position = "top")

Summary comparison

summary_gauss <- result_gauss$summary %>%
  mutate(
    prior_range = prior_upper[parameter] - prior_lower[parameter],
    post_range  = upper_90 - lower_90,
    uncert_reduction = round(100 * (1 - post_range / prior_range), 1)
  )

cat("=== Uniform Prior: Posterior Summary ===\n")
## === Uniform Prior: Posterior Summary ===
cat(sprintf("%-8s  %10s  %10s  %10s\n", "Param", "Median", "90% lo", "90% hi"))
## Param         Median      90% lo      90% hi
for (i in seq_len(nrow(summary_with_reduction))) {
  r <- summary_with_reduction[i, ]
  cat(sprintf("%-8s  %10.4f  %10.4f  %10.4f\n",
              r$parameter, r$median, r$lower_90, r$upper_90))
}
## a21           0.1057      0.0695      0.1573
## a32           0.0119      0.0015      0.0317
## input        46.9212     16.7631    106.4977
## k1            0.0323      0.0081      0.0635
## k2            0.0188      0.0060      0.0348
## k3            0.0256      0.0040      0.0518
cat("\n=== Gaussian Prior: Posterior Summary ===\n")
## 
## === Gaussian Prior: Posterior Summary ===
cat(sprintf("%-8s  %10s  %10s  %10s\n", "Param", "Median", "90% lo", "90% hi"))
## Param         Median      90% lo      90% hi
for (i in seq_len(nrow(summary_gauss))) {
  r <- summary_gauss[i, ]
  cat(sprintf("%-8s  %10.4f  %10.4f  %10.4f\n",
              r$parameter, r$median, r$lower_90, r$upper_90))
}
## a21           0.0907      0.0566      0.1285
## a32           0.0078      0.0009      0.0261
## input        50.1181     18.6064    106.6724
## k1            0.0307      0.0123      0.0553
## k2            0.0180      0.0058      0.0308
## k3            0.0115      0.0018      0.0236

What to notice:

  • For well-constrained parameters (those where the uniform posterior was already narrow relative to the prior), the Gaussian prior makes little difference. The data are strong enough to override the prior shape.
  • For weakly constrained parameters (where the uniform posterior filled much of the prior range), the Gaussian prior can noticeably tighten the posterior. The bell-shaped prior gently pulls the posterior toward the literature-informed center.
  • The Gaussian prior can also shift the posterior slightly if the literature center doesn’t perfectly align with the data – this is prior influence at work, and whether it’s helpful depends on whether the literature values are actually applicable to your site.

The practical takeaway: when data are limited (as they often are in soil science), the shape of your prior is part of your scientific argument. Gaussian priors informed by literature values are a principled way to incorporate existing knowledge. But always run a sensitivity comparison like this one to verify that your conclusions aren’t driven by the prior alone.

10. Comparing 1-Pool and 3-Pool Results

How does the 3-pool model compare to the 1-pool model? Let’s run a quick 1-pool calibration for comparison and look at bulk predictions from both models.

# Aggregate to bulk observations (same as 1-pool tutorial)
dat_1pool <- datComb %>%
  filter(year == 2009 | treatment == "Control") %>%
  mutate(C_stock_g    = C_stock * 1000,
         C_stock_sd_g = C_stock_sd * 1000) %>%
  group_by(year) %>%
  summarise(
    C14        = weighted.mean(C14, C_stock_g),
    C_14_sd    = sqrt(sum((C_stock_g * C_14_sd)^2)) / sum(C_stock_g),
    C_stock    = sum(C_stock_g),
    C_stock_sd = sqrt(sum(C_stock_sd_g^2)),
    .groups = "drop"
  )

C0_1pool <- dat_1pool$C_stock[dat_1pool$year == 2009]
F0_1pool <- dat_1pool$C14[dat_1pool$year == 2009]
years_1pool <- seq(2009, 2022)
obs_2022_1pool <- dat_1pool %>% filter(year == 2022)

likelihood_1pool <- function(pars) {
  k     <- pars[1]
  input <- pars[2]
  model <- tryCatch(
    OnepModel14(t = years_1pool, k = k, C0 = C0_1pool,
                F0_Delta14C = F0_1pool, In = input,
                inputFc = atmIn, pass = TRUE),
    error = function(e) return(NULL)
  )
  if (is.null(model)) return(-1e10)
  pred_C14    <- getF14(model)[length(years_1pool), ]
  pred_Cstock <- getC(model)[length(years_1pool), ]
  ll_C14    <- dnorm(pred_C14 - obs_2022_1pool$C14, 0,
                     obs_2022_1pool$C_14_sd, log = TRUE)
  ll_Cstock <- dnorm(pred_Cstock - obs_2022_1pool$C_stock, 0,
                     obs_2022_1pool$C_stock_sd, log = TRUE)
  return(sum(ll_C14, ll_Cstock, na.rm = TRUE))
}

result_1pool <- run_and_summarize(
  likelihood_fn = likelihood_1pool,
  prior_lower   = c(k = 0.001, input = 50),
  prior_upper   = c(k = 0.1,   input = 1000),
  par_names     = c("k", "input"),
  iterations    = 10000,
  burn_in       = 2000,
  thin          = 5
)
years_plot_cmp <- seq(2009, 2022, by = 0.5)

# 1-pool MAP prediction
map_1pool <- OnepModel14(
  t = years_plot_cmp, k = result_1pool$map["k"], C0 = C0_1pool,
  F0_Delta14C = F0_1pool, In = result_1pool$map["input"],
  inputFc = atmIn, pass = TRUE
)
pred_1pool <- data.frame(
  year = years_plot_cmp,
  C14 = as.numeric(getF14(map_1pool)),
  C_stock = as.numeric(getC(map_1pool)),
  model = "1-pool"
)

# 3-pool MAP prediction, aggregated to bulk
map_3pool_model <- ThreepSeriesModel14(
  t = years_plot_cmp,
  ks = c(k1 = result$map["k1"], k2 = result$map["k2"], k3 = result$map["k3"]),
  C0 = C0_3pool, F0_Delta14C = F0_3pool,
  In = result$map["input"], a21 = result$map["a21"], a32 = result$map["a32"],
  inputFc = atmIn, pass = TRUE
)
C_3pool   <- getC(map_3pool_model)
C14_3pool <- getF14(map_3pool_model)
# Bulk = sum of C stocks; bulk 14C = stock-weighted mean
bulk_C     <- rowSums(C_3pool)
bulk_C14   <- rowSums(C_3pool * C14_3pool) / bulk_C

pred_3pool <- data.frame(
  year = years_plot_cmp,
  C14 = as.numeric(bulk_C14),
  C_stock = as.numeric(bulk_C),
  model = "3-pool (bulk)"
)

pred_compare <- bind_rows(pred_1pool, pred_3pool)

p_cmp_c14 <- ggplot(pred_compare, aes(x = year, y = C14, color = model)) +
  geom_line(linewidth = 1.2) +
  geom_point(data = dat_1pool, aes(x = year, y = C14),
             inherit.aes = FALSE, size = 4) +
  geom_errorbar(data = dat_1pool,
                aes(x = year, ymin = C14 - C_14_sd, ymax = C14 + C_14_sd),
                inherit.aes = FALSE, width = 0.3) +
  scale_color_manual(values = c("steelblue", "coral")) +
  labs(title = expression(paste("Bulk ", Delta^14, "C: 1-Pool vs. 3-Pool")),
       x = "Year", y = expression(paste(Delta^14, "C (per mil)")),
       color = NULL) +
  theme_minimal(base_size = 14) +
  theme(legend.position = "top")

p_cmp_cstock <- ggplot(pred_compare, aes(x = year, y = C_stock / 1000, color = model)) +
  geom_line(linewidth = 1.2) +
  geom_point(data = dat_1pool, aes(x = year, y = C_stock / 1000),
             inherit.aes = FALSE, size = 4) +
  geom_errorbar(data = dat_1pool,
                aes(x = year, ymin = (C_stock - C_stock_sd) / 1000,
                    ymax = (C_stock + C_stock_sd) / 1000),
                inherit.aes = FALSE, width = 0.3) +
  scale_color_manual(values = c("steelblue", "coral")) +
  labs(title = "Bulk Carbon Stock: 1-Pool vs. 3-Pool",
       x = "Year", y = expression(paste("C stock (kg m"^-2, ")")),
       color = NULL) +
  theme_minimal(base_size = 14) +
  theme(legend.position = "top")

p_cmp_c14 / p_cmp_cstock

What does the 3-pool model tell us that the 1-pool model can’t?

  • Layer-specific turnover times: The 3-pool model estimates that surface organic carbon cycles in ~45 years, while mineral carbon takes ~81 years. The 1-pool model gives a single bulk estimate that averages over these very different timescales.

  • Carbon transfer through the profile: The transfer rates \(a_{21}\) and \(a_{32}\) quantify how fast carbon moves downward from surface to deep to mineral soil. This is invisible to the 1-pool model.

  • Layer-specific predictions: The 3-pool model predicts how each individual layer will evolve, not just the bulk. This matters for understanding vulnerability – surface organic carbon may respond very differently to warming than deep mineral carbon.

The cost of complexity: More parameters mean wider posteriors (less certainty per parameter), stronger correlations between parameters, and harder convergence. The 6-parameter posterior is much harder for MCMC to explore than the 2-parameter posterior. This is the classic bias-variance tradeoff in modeling – more complex models can capture more reality, but they’re harder to constrain from limited data.

11. Wrapping Up

Key takeaways

  1. Multi-pool models capture real structure that single-pool models miss. The three soil layers have genuinely different turnover times, and a 3-pool model can estimate these separately.

  2. More parameters = harder inference. With 6 parameters and 6 data points, parameter correlations become important. Always check the posterior correlation matrix – it tells you which parameters the data can and can’t independently constrain.

  3. Prior shape matters for weakly constrained parameters. We saw in Section 9 that switching from uniform to Gaussian priors can noticeably affect posteriors for parameters with weak data constraint. When data are limited, the prior is part of your scientific argument – choose it deliberately and always run a sensitivity analysis.

  4. Trace plot diagnosis is even more critical with more parameters. Some parameters (especially for slow pools) may mix poorly because the data provide weak constraints. This is a feature, not a bug – it tells you where the data are uninformative.

  5. The uncertainty scaling factor (0.4) on the likelihood matters a lot. It accounts for the gap between measurement uncertainty and model structural error. Experiment with it to understand its effect.

  6. Model comparison is essential. Comparing 1-pool and 3-pool predictions side by side (Section 10) reveals when the extra complexity is justified and when a simpler model suffices.

Suggestions for further exploration

  • Try feedback topology: Instead of a pure series model (carbon only flows downward), allow some upward transfer. How does this change the posteriors?

  • Try the warming treatment: The datComb table also contains a “Warming” treatment. What happens to decomposition rates under experimental warming?

  • Try log-normal priors for k: Decomposition rates span orders of magnitude, making log-normal priors a natural choice. Use createPrior() with a custom density and sampler function.

  • Increase iterations: Try 100,000 iterations. Do the posteriors change? If so, 20,000 wasn’t enough for full convergence.

sessionInfo()
## R version 4.4.0 (2024-04-24)
## Platform: aarch64-apple-darwin20
## Running under: macOS 26.0.1
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.0
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: America/Phoenix
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] patchwork_1.3.2     tidyr_1.3.1         dplyr_1.1.4        
## [4] ggplot2_4.0.0       BayesianTools_0.1.8 SoilR_1.2.107      
## [7] deSolve_1.40       
## 
## loaded via a namespace (and not attached):
##  [1] gtable_0.3.6         xfun_0.52            bslib_0.7.0         
##  [4] sets_1.0-25          lattice_0.22-6       vctrs_0.6.5         
##  [7] tools_4.4.0          Rdpack_2.6.6         generics_0.1.3      
## [10] parallel_4.4.0       tibble_3.2.1         fansi_1.0.6         
## [13] highr_0.10           pkgconfig_2.0.3      Matrix_1.7-0        
## [16] DHARMa_0.4.7         RColorBrewer_1.1-3   S7_0.2.0            
## [19] assertthat_0.2.1     lifecycle_1.0.4      compiler_4.4.0      
## [22] farver_2.1.2         stringr_1.5.1        Brobdingnag_1.2-9   
## [25] htmltools_0.5.8.1    sass_0.4.9           yaml_2.3.8          
## [28] pillar_1.9.0         nloptr_2.2.1         jquerylib_0.1.4     
## [31] MASS_7.3-60.2        cachem_1.1.0         reformulas_0.4.4    
## [34] bridgesampling_1.2-1 boot_1.3-30          nlme_3.1-164        
## [37] tidyselect_1.2.1     digest_0.6.35        mvtnorm_1.2-4       
## [40] stringi_1.8.4        purrr_1.1.0          labeling_0.4.3      
## [43] splines_4.4.0        fastmap_1.2.0        grid_4.4.0          
## [46] expm_1.0-0           cli_3.6.2            magrittr_2.0.3      
## [49] survival_3.5-8       utf8_1.2.4           withr_3.0.0         
## [52] scales_1.4.0         rmarkdown_2.29       igraph_2.1.4        
## [55] lme4_1.1-38          msm_1.8.2            coda_0.19-4.1       
## [58] evaluate_0.23        knitr_1.46           rbibutils_2.4.1     
## [61] rlang_1.1.3          Rcpp_1.0.12          glue_1.7.0          
## [64] minqa_1.2.8          jsonlite_1.8.8       R6_2.5.1